Human-caused wolf mortality persists for years after discontinuation of hunting

By the mid-twentieth century, wolves were nearly extinct in the lower 48 states, with a small number surviving in northern Minnesota. After wolves were placed on the endangered species list in 1973, the northern Minnesota wolf population increased and stabilized by the early 2000s. A wolf trophy hunt was introduced in 2012–2014 and then halted by a court order in December 2014. The Minnesota Department of Natural Resources collected wolf radiotelemetry data for the years 2004–2019. Statistical analysis showed that wolf mortality remained close to constant from 2004 until the initiation of the hunt, and that mortality doubled with the initiation of the first hunting and trapping season in 2012, remaining at a nearly constant elevated level through 2019. Notably, average annual wolf mortality increased from 21.7% before wolf hunting seasons (10.0% by human causes and 11.7% natural causes) to 43.4% (35.8% by human causes and 7.6% natural causes). The fine-grained statistical trend implies that human-caused mortality increased sharply during the hunting seasons, while natural mortality initially dropped. After the hunt’s discontinuation, human-caused mortality remained higher than prior to the hunting seasons throughout the five years of the available after-hunt radiotelemetry data.

The R-days column title stands for radiotelemetry days. Fractional wolf death counts are the result of six wolf deaths with unknown causes, imputed proportionately to the known causes of mortality for each period.
Supplementary  Table 2 has 0 wolf deaths (no variance) and quite a low number of radio-days, with one wolf death producing an outlier in the year 2019. Although the outliers were handled well by the variance weighting, the zero variance required an additional estimation of possible variance, perhaps from nearby valid entries. Additionally, with high-granularity radiotelemetry (as used in MN DNR data collection), the greater variance in years with higher wolf mortality may create its own bias. A second weighting approach involves using radio-days: where V is variance, p and (1 − p) are the probability fractions, D is the number of radio-days (measurements), and K is the count of wolf deaths in the radio-days period. With a large D and small K, (D − K)/D ≈ 1, and K/D 2 ∼ 1/D 2 . The final step shows the use of the unbiased number of recorded deaths weighted by radio-days; thus, when weighting can be performed as 1/V , it must be performed alternatively as D 2 . The advantage of this second weighting approach is its consistency, and it does not require inventing a special treatment for the years with no recorded wolf mortality. The following comparisons show that this weighting method also allows a more precise match of the regression to the summary analysis in Table 1. Finally, the R models for these two weights are implemented as follows: and where wFit1 and wFit2 are the corresponding weighted fit scenarios for the linear regression, Y are the typical daily hazard values or survival, X is the year range, D is the number of radio-days (measurements) per year, and df is an R data frame with a loaded CSV input (see also Sheather [1]). On the following three pages, the tables are grouped in comparison pairs for weighting using inverse variance and D 2 (see Supplementary  Table 3-Supplementary Table 4, Supplementary Table 5-Supplementary  Table 6, and Supplementary Table 7-Supplementary Table 8). For completeness, these tables show the parameters of the linear regression with years offset to the first day of each month of the year. The most important columns are stdError (the standard error of the regression), r 2 (the goodness-of-fit), and mean survival, which is typically calculated as: where SurvM is the mean survival over the subsequent number of years T in the regression interval Y min to Y max . As shown below, stdErr is comparatively lowest for years offset to start as of November 1st in all tables except Supplementary Table 4, where it is a close second. The regression parameter r 2 is a goodness-of-fit of the trend. In Supplementary Table 3-Supplementary  Table 6 r 2 is low, representing a near-constant level of survival and mortality. For the entire dataset duration the best regression was also observed with a November 1st time frame offset. The second equally important consideration for the weighting method is how well the mean survivals in these tables match the multiyear summaries in Table 1. Here, the summaries from sections in Table 1  Consequently, the regression analysis when using inverse variance weighting resulted in nearly identical but slightly underestimated mortality trend values. Thus, we used the weighting by radio-days squared in our reported analysis.

Supplementary Note 3. Pinpointing the trend change timing via regression discontinuity analysis
Regression discontinuity analysis was performed to determine whether there was a cutoff point where the trend may have become discontinuous [2], as could occur before and after a specific date or event, with trends having the potential to significantly change in their magnitude and slope. This was performed using the R script WolfProc.R as a varying slope regression discontinuity analysis [3][4][5].
Supplementary The R script iterated through each month of five consecutive years and tested for the existence of discontinuity cutoff points, see Supplementary  Table 9. The only sharp trend discontinuity cutoff [2]  The January row in Supplementary Table 9 shows that statistics aligned with standard calendar years starting on January 1st were particularly unsuitable for noticing sharp discontinuity. Only part of the mortality from the 2012 to 2013 hunting season occurred in November and December 2012, thus producing a smooth transition when counting statistics in calendar years. Supplementary Fig. 1 graphically illustrates the numerical p-values in Supplementary Table 9. Testing for the cutoff points on November 1st of years 2011 and 2013 showed a low likelihood of such cut point, with p-values of 0.813 and 0.684, respectively, leading to the appearance of Supplementary  Fig. 1a and Supplementary Fig. 1c. On the other hand, the p-value was 0.039 for November 2012 and corresponds to a sharp breaking of the trend in Supplementary Fig. 1b. There is always an uncertainty inherent to radiotelemetry data when accounting for the fate of the censored wolves and the exact dating of the deaths that may lead to mismeasured wolf mortality [6]. Our analysis used the end date in the MN DNR data for the established wolf death date. There may be other choices of adjusting the uncertainty of the death date, as was done by Chakrabarti et al. [7], however we found that for this data set such adjustments have minimal effect on the survival/mortality results, and any adjustment method without the precise timing knowledge remains arbitrary. A validation using the last seen alive date as the mortality date was performed, demonstrating that even with such extreme edge case adjustment the overall mortality would change by only a fraction of a percent (from 34.7% to 35.6%), and confirming a similar trend discontinuity in October-November 2012 (see in Supplementary Table 10). Thus, the end date in MN DNR data set was chosen for consistency of treatment with the censored wolves.
A verification for the trend discontinuity within the during-and-after hunt data intervals is performed in Supplementary Note 4.   Table 12 presents the regression analysis performed while discounting 6 out of 59 wolf mortality events due to unknown causes, rather than imputing them as fractional shares of mortalities. Imputing the unknown cause mortality allows us to account for all wolf deaths in the data rather than underestimating wolf mortality, with no noticeable bias in the regression outcome.
Supplementary See the comparisons between Supplementary Table 12 and corresponding  Table 1 (where imputation was performed), which show that the patterns appear qualitatively the same between the two representations, with 10% lower values across all numbers where the unknown deaths were omitted, thereby showing the benefit of this imputation.